ORIGINAL 


ARTICLE 


Asian Herpetological Research 2021, 12(3): 295-307 
DOI: 10.16373/j.cnki.ahr.200129 


Complete Mitochondrial Genome of the Steppe Ribbon Racer (Psammophis 
lineolatus): The First Representative from the Snake Family Psammophiidae 


and its Phylogenetic Implications 


Minli CHEN", Jinlong LIU’, Jun LI’, Na WU" and Xianguang GUO” 


"Chengdu Institute of Biology, Chinese Academy of Sciences, Chengdu 610041, Sichuan, China 


“University of Chinese Academy of Sciences, Beijing 100049, China 


* College of Life Science and Technology, Xinjiang University, Urumqi 830046, Xinjiang, China 


Abstract Comparisons of vertebrate mitochondrial 
genomes (mitogenomes) may yield significant insights 
into the evolution of organisms and genomes. However, 
no complete mitogenome from the snake family 
Psammophiidae has been reported. In this study, we 
sequenced and annotated the complete mitogenome of 
Psammophis lineolatus, representing the first mitogenome 
of Psammophiidae. The total length is 17 166 bp, 
consisting of 13 protein-coding genes (PCGs), 22 transfer 
RNAs (tRNAs), two ribosomal RNAs (12S rRNA and 
16S rRNA), and duplicate control regions (CRs). This 
gene arrangement belongs to the Type III pattern, 
which is a widely shared gene order in Alethinophidian 
snakes. All tRNAs exhibit cloverleaf structures with 
the exception of tRNA-Ser*”* and tRNA-Pro, which 
lack a dihydrouridine (DHU) arm/stem and TYC loop, 
respectively. The 13 PCGs include five start codons (ATG, 
GTG, ATA, ATT, and ATC), two complete stop codons 
(TAA and AGG), and two incomplete stop codon (T-- 
and TA-). In addition, the Ka/Ks ratios indicate that 
all PCGs had undergone a strong purifying selection. 
Four types of CR domains rearrangement occurred 
among eight species of Elapoidea. The phylogenetic 
reconstructions with both Bayesian inference and 
maximum likelihood methods support the placement 
of Psammophiidae in the Elapoidea superfamily, 
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with Homalopsidae being the sister taxon to Elapoidea 
and Colubroidea. However, the sister taxon of 
Psammophiidae is unclear due to the availability of 
Elapoidea mitogenomes being limited to the family 
Elapidae. More mitogenomes from different taxonomic 
groups in Elapoidea are needed to better understand the 
phylogenetic relationships within Elapoidea. 


Keywords mitochondrial genome, phylogenetic relationships, 
Psammophis lineolatus, Pammophiidae, purifying selection 


1. Introduction 


During the last four decades, sequences of the mitochondrial 
genomes (mitogenomes) of thousands of Metazoan species have 
been produced and are freely available in GenBank (https:// 
www.ncbi.nlm.nih.gov/genbank/). Approximately two thirds 
of these sequences are from vertebrates. These sequences have 
been successfully applied to phylogenetic analyses at both 
shallow and deep levels (e.g., Bernt et al., 2013a; Knaus et al, 
2011), as well as for analyses of mitogenome organization and 
evolution (eg, Bernt et al, 2013b; Boore, 1999; Mauro et al, 2006). 
As the number of mitogenome sequences continues to increase, 
it has become clear that, rearrangements of several transfer 
RNA genes (tRNAs) are quite common (eg, Mauro et al, 2006; 
Pereira, 2000; Xia et al, 2016), even though gene arrangement is 
quite stable across many vertebrates. Overall, most vertebrate 
mitogenomes contain the same 37 genes (two for rRNAs, 13 
for proteins and 22 for tRNAs), but their order is variable 


among - and to a lesser extent, within — major groups (Bernt 
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et al, 2013; Boore, 1999; Pereira, 2000). Accordingly, comparisons 
of vertebrate mitochondrial genomes may lead to significant 
insights into the evolution both of organisms and of genomes. 

Snakes, which are members of the suborder Serpentes in 
Squamata, are an example of a group having a dynamic 
mitogenome structure. Currently, there are over 3700 
recognized snake species, comprising two major clades: 
Scolecophidia (blind snakes) and Alethinophidia (Uetz et al., 
2020). Scolecophidia are small fossorial snakes that feed mainly 
on ants and termites. Alethinophidia includes Henophidia 
(primitive snakes) and Caenophidia (advanced snakes; including 
~80% of extant snakes) (Vidal and Hedges, 2002a). Some unique 
characteristics have been reported in the mitogenomes of 
snakes, including duplicate control regions (CRs), gene order 
rearrangements, shorter tRNA genes, and other shortened 
genes (Chen and Zhao, 2009; Douglas and Gower, 2010; Jiang 
et al, 2007). To date, 11 types of rearrangement patterns have 
been detected in snake mitogenomes (Type I, Il, III, TILA, III-B, 
IN-B1, TIC, I-D, TEE, MEF, III-G) (Qian et al., 2018). Among 
these, Type III is characterized by duplication of the CR and 
translocation of the tRNA-Leu gene; it is also the most common 
gene arrangement, extensively distributed in Alethinophidia. 
These arrangements are thought to involve three main 
processes: gene loss, translocation and duplication (Qian et al, 
2018). Although a few convergent gene rearrangements have 
been observed in vertebrate mtDNAs, many derived genomic 
rearrangements are clearly useful as phylogenetic makers 
(thus, synapomorphic characters) for the identification of 
monophyletic groups (Boore and Brown, 1998; Okajima and 
Kumazawa, 2010; Xia et al, 2014). Therefore, mitogenome data 
will contribute to our phylogenetic understanding of snakes, 
as well as the evolutionary implications of their mitogenomic 
rearrangements. 

As of October 1, 2020, there were 188 complete mitogenomes 
of snakes present in the NCBI RefSeq database (Pruitt et al., 
2007), covering 124 snake species. However, no complete 
mitogenome from the family Psammophiidae (formerly the 
subfamily Psammophiinae; Bourgeois, 1968) has ever been 
reported. This family comprises a natural group occurring 
throughout Africa, the Middle East, Madagascar, southern 
Europe and south-central Asia (Dowling and Duellman, 1978; 
Kelly et al., 2003; Kelly et al, 2009), and currently consists of 
eight genera and 55 extant species (Uetz et al, 2020). Many 
studies of Psammophiidae have focused on systematics 
based on morphological features (Bourgeois, 1968; Dowling 
and Duellman, 1978; Zaher, 1999). Although there are some 
molecular studies, most are based on limited mitochondrial 
DNA fragments and/or nuclear genes (Gravlun, 2001; Kelly et 
al, 2008; Kelly et al., 2009; Lawson et al, 2005; Nagy et al, 2003; 
Vidal and Hedges, 2002b). 
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In this study, we produced for the first time the complete 
mitogenome sequence of the Steppe Ribbon Racer, Psammophis 
lineolatus (Brandt, 1838), a representative of the family 
Psammophiidae. As noted by Qian et al. (2018), the Type 
III mitogenome rearrangement is the most prevalent and 
widely occurring in Elapidae. Given that Psammophiidae is 
closely related to Elapidae among the major snake families 
- both belonging to the Elapoidea superfamily (Figueroa et 
al, 2016; Kelly et al, 2009; Zaher et al, 2019) - we predict that 
the mitogenome arrangement of P. lineolatus would belong to 
the Type III pattern. In addition, we used the mitogenome to 
infer the phylogenetic position of Psammophiidae in the major 
lineages of the advanced snakes. 


2. Materials and Methods 


2.1. Sample collection and DNA extraction A P. lineolatus 
specimen was captured by hand in Buerjin County (47.559316°N, 
87.056509°E), Xinjiang Uygur Autonomous Region, China in 
July 2019. A piece of liver was dissected from the snake, which 
has been euthanized with an overdose of sodium pentobarbital 
delivered by intraperitoneal injection. Both the liver sample 
and voucher specimen (Field number GXG311) were fixed in 
95% ethanol, and deposited at the herpetological collection, 
Chengdu Institute of Biology, Chinese Academy of Sciences. 
Total genomic DNA was extracted using the TIANamp 
Genomic DNA Kit (TIANGEN, Beijing, China) following the 
manufacturers instructions, and stored at -20°C. The amount 
and quality of the extracted DNA were assessed through 
electrophoresis in a 0.8% agarose gel with GoldView (Solarbio, 
Beijing, China). 


2.2. PCR amplification and sequencing Twenty pairs of 
primers (Table SI) were designed to amplify the complete 
mitogenome, based on published mitogenomes from Elapidae 
(GenBank accession number: EU579523, GU045453). The 
polymerase chain reaction (PCR) amplification was carried out 
in volumes of 25 uL consisting of 2 uL template DNA (~ 40 ng/ 
uL), 25 yL dNTPs (02 mmol/L each), 15 yL MgSO, (1.5 mmol/ 
L), 2.5 uL buffer of KOD-plus neo (TOYOBO, Shanghai, China), 
05 pL KOD-plus neo (1.0 U/L), 0.75 pL each primer (0.3 umol/ 
L), and 145 pL ultrapure water. Thermal cycling was performed 
with an initial denaturation at 94°C for 2 min, followed by 32- 
35 cycles of 98°C for 10 s, 45-60°C (depending on the primers) 
for 30 s, and 68°C for 60-90 s; the reaction concluded with a 
final extension of 7 min at 68°C, then stored at 4°C. The PCR 
products were assessed using 1% agarose gel electrophoresis, 
purified, and sequenced with the PCR primers and internal 
primers generated by primer walking. All fragments were 
sequenced on an ABI 3730 automated sequencer (Applied 
Biosystems, Inc.) at TsingKe Biotech (Beijing, China). 
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2.3. Genome assembly, annotation and analysis 
Mitochondrial DNA fragments were assembled to create the 
complete mitochondrial genome using the SeqMan II program 
included in the Lasergene software package (DNAStar Inc., 
Madison, USA). The complete mitogenome was submitted to 
GenBank using the software of Sequin v15.10 (http://www. 
ncbi.nlm.nih.gov/Sequin/). The boundaries of the protein- 
coding genes (PCGs) and ribosomal RNA genes (rRNAs) were 
identified using NCBI-BLAST (http://blast.ncbi.nlm.nih.gov) 
and MITOS Web Server (http://mitos2.bioinf.uni-leipzig.de/ 
index.py) (Bernt et al., 2013c). The tRNAs and their potential 
cloverleaf structures were predicted using tR NAscan-SE Search 
Server v1.21 (http://lowelab.ucsc.edu/tRNAscan-SE/) (Lowe and 
Chan, 2016), and verified by aligning them with other snakes. 
The size of the CR was confirmed by the boundaries of tRNAs 
and by sequence comparison with previously reported snake 
mitogenomes. A circular mitogenomic map was generated 
using OGDRAW v1.3.1 (Greiner et al, 2019). Base composition 
and the relative synonymous codon usage (RSCU) values 
were calculated using MEGA v7.0 (Kumar et al., 2016). Strand 
asymmetry was calculated using the formulae: AT-skew = (A 
- DAA + T); GC-skew = (G - C)/(G + C) (Perna and Kocher, 
1995). 

Ka/Ks is the ratio of the number of nonsynonymous 
substitutions per nonsynonymous site (Ka) to the number of 
synonymous substitutions per synonymous site (Ks) (Hurst, 
2002). This ratio is considered to be a good indicator of selective 
pressure at the sequence level (Yang and Bielawski, 2000) and 
has been used to identify PCGs under positive and purifying 
selection in a wide range of organisms (Hurst, 2009). In 
general, positive selection yields a Ka/Ks ratio greater than 1. 
Conversely, if purifying selection occurs, the Ka/Ks ratio is less 
than 1, which would lead to a decrease in diversity at the amino 
acid level. Therefore, to determine the direction of natural 
selection (positive or purifying) of the PCGs at the molecular 
level, we calculated the Ka/Ks ratio of each PCG between P. 
lineolatus and three Elapidae snakes (Bungarus multicinctus, Naja 
naja and Ophiophagus Hannah) using DnaSP v6.12.03 (Rozas et 
al, 2017). 


2.4. Phylogenetic reconstruction Phylogenetic analyses 
were conducted using the mitochondrial genome (excluding 
the CR) to infer the placement of Psammophiidae within 
the advanced snakes (Table S2). Complete gene sequences 
from representatives of the major snake families were chosen 
following a recent study (Qian et al., 2018), and were retrieved 
from GenBank. Each gene was extracted from the whole 
mitogenome, then separately aligned in batches with MAFFT 
v7.263 (Katoh and Standley, 2013) using ‘auto’ strategy and 
normal alignment mode. In total, 37 mitochondrial genes were 
combined from 71 taxa using SeaView v5.0.1 (Gouy et al, 2010), 
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which formed a dataset of 15 838 bp in length. Considering 
that Henophidia is closely related to Caenophidia (Streicher 
and Wiens, 2016; Yan et al, 2008), six Henophidia snake species, 
representing six families, were chosen as outgroup taxa (Table 
S2). Among these, Tropidophiidae and Aniliidae were used to 
root the trees based on two recent publications Qian et al, 2018; 
Streicher and Wiens, 2016). 

Prior to the phylogenetic analyses, optimal combinations of 
partitioning schemes and nucleotide substitution models were 
determined with PartitionFinder v2.11 (Lanfear et al., 2017), 
using the “greedy” algorithm and the Akaike information 
criterion. The branch lengths of alternative partitions were 
linked and the “mrbayes” model was used to estimate the best 
substitution model for IQ-TREE and MrBayes. The best-fit 
substitution models and partitioning schemes for each gene are 
given in Table S3. 

Phylogenetic analyses were performed with Bayesian 
inference (BI) and maximum likelihood (ML) methods. Bayesian 
phylogenies were inferred using MrBayes v3.2.6 (Ronquist et al, 
2012) under partitioned models. Two independent runs were 
conducted with four Markov Chain Monte Carlo (MCMC) for 
10,000,000 generations with parameters and topologies sampled 
every 1000 generations. When the average standard deviation 
of split frequencies reached a value of less than 0.01, the initial 
25% trees were discarded as burn-in and the remaining trees 
were used to calculate Bayesian posterior proba bilities (PP). 
Maximum likelihood analyses were conducted in IQ-TREE 
v16.12 (Nguyen et al., 2015) under the models selected for each 
identified partition. We assessed nodal support using 5000 
bootstrap pseudoreplicates via the ultrafast-bootstrap (UFBoot) 
(Minh et al., 2013) method. Nodes with UFBoot 2 95 were 
considered to be well-supported (Minh et al, 2013; Wilcox et al, 
2002). Finally, the trees were visualized and edited in FigTree 
v1.42 (http://tree.bio.ed.acuk/software/figtree) and PowerPoint. 


3. Results and Discussion 


3.1. Genome structure, organization and composition The 
complete sequence of the mitogenome of P. lineolatus (accession 
number MT991050) is a traditional circular, double- strain 
DNA molecule (Figure 1). It contains 37 genes with 13 PCGs, 
two rRNA genes, and 22 tRNA genes but two control regions 
(CR I, CR II) (Table 1). However, the tRNA-Leu""" (Leu2) gene 
was translocated; in most vertebrates, it is between 16S rRNA 
and NDI, but in snakes it is between tRNA-Ile and tRNA-GIn 
(Figure 1). The gene order is in agreement with the Type III 
pattern designated by Qian et al. (2018). The majority of the 37 
genes were encoded by the heavy-strand, except the ND6 gene 
and eight tRNAs (tRNA-GIn, Ala, Asn, Cys, Tyr, Ser2, Glu, Pro). 
The intergenic spacers were 33 bp, ranging from 1 to 9 bp in 
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Figure 1 Mitochondrial map of Psammophis lineolatus. Genes encoded on the heavy or light strand are indicated on the outside or 
inside of the circular mitogenome map, respectively, showing the direction of transcription. The tRNAs are denoted by color and labeled 


according to the three-letter amino acid codes. 


size; the longest (9 bp) was located between ND6 and tRNA- 
Glu. Overlapping nucleotides were found between six pairs of 
genes, with the length varying from 1 to 25 bp; the largest was 
located between tRNA-Val and 16S rRNA. 

The whole mitogenome of P. lineolatus is 17 166 bp in 
length (Table 1), with A 32%, G 15.5%, T 23.596, C 29% (Table 
2). The strand bias of nucleotide composition of metazoans 
mtiogenomes are usually measured by AT-skews and GC- 
skews (Perna and Kocher, 1995). The AT-skew and GC-skew 
of the mitogenome are 0.154 and -0.305, respectively (Table 2), 
indicating As and Cs are more abundant than Ts and Gs. The 
nucleotide composition of P. lineolatus is obviously inclined to 
A/T. 


3.2. Protein-coding genes and codon usage The total 
length of the PCGs in mitogenome of P. lineolatus was 11 255 
bp, which consisted of seven NADH dehydrogenases (NDI- 
ND6 and ND4L), three cytochrome c oxidases (COX1-COX3), 
two ATPases (ATP6 and ATP8) and one cytochrome b (Cyt 
b) (Table 1). Herein, five-type starting codons for PCGs were 
adopted: most start with ATG (COX2, ATP8, ATP6, COX3, 
NDAL, ND4, ND5, ND6, Cytb), while NDI with ATA, ND2 
with ATT, COX1 with GTG and ND3 with ATC. The majority 
of the PCGs terminated with TAA, while COX1 used AGG. 
Five (ND2, COX2, COX3, ND3 and Cytb) terminated with 
incomplete T, presumably transformed into complete stop 
codons through post-transcriptional polyadenylation (Ojala et 
al, 1981). 
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Table 1 Complete mitochondrial genome of Psammophis lineolatus and traits of codon as well as anticodon. 


Position Intergenic Codon . 

Gene/element From To ducleotde Length (bp) San Stop Anticodon Stand 
tRNA-Phe 1 66 0 66 GAA H 
12S rRNA 67 991 0 925 H 
tRNA-Val 989 1053 ES 65 TAG H 
16S rRNA 1029 2537 -25 1509 H 
ND1 2546 3494 8 949 ATA Tas H 
tRNA-Ile 3495 3559 0 65 GAT H 
CRII 3560 4571 0 1012 H 
tRNA-Leu ""(Leu2) 4572 4644 0 73 TAA H 
tRNA-GIn 4645 4715 0 71 TIG iL 
tRNA-Met 4717 4779 1 63 CAT H 
ND2 4780 5806 0 1027 ATT Wee H 
tRNA-Trp 5807 5870 0 64 TCA H 
tRNA-Ala 5873 5938 3 66 NGC L 
tRNA-Asn 5939 6009 0 71 GTT L 
O, 6012 6049 2 38 

tRNA-Cys 6047 6107 3 61 GCA L 
tRNA-Tyr 6109 6169 1 61 GTA iL, 
COXI 6171 7772 1 1602 GTG AGG H 
tRNA-Ser ^ "(Ser2) 7763 7829 0 67 TGA iL, 
tRNA-Asp 7830 7895 0 66 GTC H 
COX2 7898 8582 2 685 ATG Tee H 
tRNA-Lys 8583 8649 0 66 TIT H 
ATP8 8653 8814 3 162 ATG TAA H 
ATP6 8805 9484 -10 680 ATG TA- H 
COX3 9485 10 268 0 784 ATG Te H 
tRNA-Gly 10 269 10 330 0 62 TCC H 
ND3 10 331 10 676 0 346 ATC Tee H 
tRNA-Arg 10 677 10 741 0 65 TCG H 
ND4L 10 742 11032 0 291 ATG TAA H 
ND4 11033 12 370 0 1338 ATG TAA H 
tRNA-His 12371 12 434 0 64 GTG H 
tRNA-Ser ?"(Ser1) 12 435 12 491 0 57 GCT H 
tRNA-Leu^ (Leu) 12 489 12 559 -3 Fil TAG H 
ND5 12 561 14 324 1 1764 ATG TAA H 
ND6 14 320 14 829 =; 510 ATG AGG iL 
tRNA-Glu 14 839 14 901 9 63 TTC E 
Cytb 14902 16 018 0 1117 ATG Te H 
tRNA-Thr 16 016 16 081 -3 66 TGT H 
tRNA-Pro 16 082 16 144 0 63 TGG L 
CRI 16 145 17 166 0 1022 H 


Thirteen PCGs contained 3 589 codons without termination 
codons. The most frequently used amino acids are Leul (7.496), 
Ser2 (6.396), and Gly (4.196), whereas the least common amino 
acids are Glu (1.8%), Cys (1.7%), Trp (0.18%) (Table 3). The RSCU 
refers to the relative frequency of each codon to code a specific 
amino acid, which could measure the bias of codon usage. The 
RSCU for the mitochondrial PCGs in P. lineolatus exhibited an 
over-usage of A and T at the third codon positions. Among 
them, AAA-Lys (1.46), CCU-Pro (1.43), ACU-Thr (1.39) and 
UCA-Ser2 (139) are the most frequently used codons. The least 


frequent codons are ACG-Thr (0.36), CCG-Pro (0.39) and GCG- 
Ala (0.49) (Table 3). The AT content of the 12 PCGs (excluding 
ND6 in the light-strand) is 54.9%; the AT-skew and GC-skew 
are 0.15 and -0.336, respectively (Table 2), which is also a similar 
trend in other snakes (Debuy et al., 2012; Wang et al, 2012). 


3.3. Transfer and ribosomal RNAs Twenty-two tRNAs were 
obtained from the complete mitochondrial of P. lineolatus. The 
tRNAs ranged from 57 to 73 bp, and the total length was 1436 
bp. There was a positive AT-skew, except in the light-strand 
tRNA (Table 1; Table 2). The skew analysis of tRNA indicated 
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Table 2 Composition and skewness of P. lineolatus mitogenome. 

Regions Strand T(U) C A G AT (%) GC (%) AT-skew GC-skew 
Whole genome H 295 29 32 15.5 559 44.5 0.154 -0.305 
PCGs H 23.3 30.1 31.6 15 54.9 45.1 0.15 -0.336 
tRNAs H 23.6 252 319 19.3 5515 44.5 0.15 -0.133 
tRNAs L 30.2 16.6 24.3 28.9 54.5 45.5 -0.109 0.269 
ATP6 H 22 30.9 31.8 15.3 53.9 46.2 0.18 —0.338 
ATP8 H 16.7 31.5 42.6 9.3 59.3 40.8 0.438 -0.545 
COXI H 255 28.1 2777 18.8 59» 46.9 0.041 —0.198 
COX2 H 22.5 30.1 27 20.4 49.5 50.5 0.091 -0.191 
COX3 H 23 30.4 29 T2 52 48.1 0.115 -0.263 
Cytb H 25.5 30.3 30.3 13.9 55.8 44.2 0.085 -0.372 
NDI H 23.8 30.1 32.5 13.6 56.3 43.7 0.154 -0.378 
ND2 H 20.5 31.6 36.2 11:7: 56.7 43.3 0.277 0.461 
ND3 H 24.3 30.9 30.6 14.2 54.9 45.1 0.116 -0.372 
ND4 H 23.7 30.6 31.8 13.9 55.5 44.5 0.147 -0.375 
ND4L H 26.1 292) 31.6 13.1 5727 42.3 0.095 -0.382 
ND5 H 22.2 30 34.5 13,2 56.7 43.2 0.217 -0.389 
ND6 IL, 47.3 8 122 32.5 595 40.5 -0.591 0.604 
16S rRNA H 19.9 24.3 374 18.1 57.6 42.4 0.308 -0.146 
12S rRNA H Tgog! 26.1 39.2 19 54.9 45.1 0.354 -0.156 
CRI H 31.9 29.2 25 14 56.9 43.2 -0.121 -0.352 
CRII H SDE DONI 24.8 14.4 56.5 43.5 -0.122 -0.338 


the use of As and Cs in the heavy-strand, while there was an 
obvious bias toward the use of Ts and Gs in the light-strand 
(Table 2). All of the tRNAs were able to fold into a typical 
cloverleaf structure, with the exception of tRNA-Ser^^' and 
tRNA-Pro, which lacked a dihydrouridine (DHU) arm/stem 
and TYC loop, respectively (Figure 2). Lacking the DHU arm/ 
stem is likely related to a compensation mechanism for tRNA 
transport function — a phenomenon that has been proposed 
in lizards (Macey et al., 1997; Yu et al, 2018) and toads (Cai et al, 
2020). The T arm contraction of snake mitochondrial tRNA 
has been reported (Kumazawa et al., 1996), and is also notable 
in P. lineolatus mitochondrial tRNA genes, most of which are 
less than 5 bp (tRNA-Ala, Asp, Glu, His, Lys, Thr) or have only 
2-3 bp (tRNA-Gly, Met, Phe). These results suggest that the 
simplification of mitochondrial tRNA is subject to evolutionary 
pressure (Kumazawa et al., 1996). 

Since there is no reported complete mitogenome in 
Psammophiidae, we identified the boundaries of rRNA genes 
based on their relative location in the mitogenome. The 12S 
rRNA and 16S rRNA genes are 925 bp and 1509 bp in length, 
respectively, and generally separated by tRNA-Val (Table 1). 
The AT content of 12S rRNA is 54.9%, and 57.6% in 16S rRNA 
of P. lineolatus (Table 2). The AT-skew of 12S rRNA and 16S 
rRNA are positive and greater than the GC-skew, suggesting 
that there are more As and Cs than Ts and Gs (Table 2). 


3.4. Control region Duplicate CRs were found in P. lineolatus 
mitogenome. The first CR (CR I, 1022 bp) was located between 
tRNA-Pro and tRNA-Phe. To date, this is the most typical 


position reported for the CR of vertebrate mtDNA. The 
second CR (CR II; 1012 bp) was located between tRNA-Ile and 
tRNA-Leu""" (Leu2) (Figure 1; Table 1). Both CRs were highly 
conserved in sequence similarity, which could be explained by 
concerted evolutionary mechanism (Kumazawa et al, 1998). All 
of their skew values were negative, indicating more Ts and Cs 
than As and Gs (Table 2). 

Three conservative functional domains were observed in 
both CRs: C-rich, termina tion-associated sequences (TASs) and 
conserved sequence blocks (CSBs). The C-rich sequence may 
facilitate formation of the D-loop by decelerating the extension 
of heavy-strand synthesis (Kumazawa et al., 1998). The CSBs 
play an important role in the CR of all organisms, which are 
associated with the start sites for DNA synthesis. It has been 
thought that TAS elements are associated with stop sites for 
DNA synthesis, and have been found near the 3' terminus of the 
D-loop (Doda et al., 1981). Generally, each control region could 
find a CSB that presents a CSB-1- CSB-2- CSB-3 arrangement. 
On the other hand, one or more TASs can appear, and can 
include forward or reverse complementary sequences. In this 
study, four types of CR domains arrangements were observed 
in the Elapoidea species (Figure 3; Table S4), and each conserved 
domain appears only once, rather than, repeatedly. Type A 
occurred in P. lineolatus; Types B and C were assigned to Naja 
spp. and O. hannah, respectively; and Type D was assigned to 
Bungarus spp. and Micrurus fulvius. These four arrangement 
types were mostly caused by the change in position of TAS2 
and the presence or absence of CSB-2. However, CSB-2 
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Table 3 Genetic code, amino acid composition and relative synonymous codon usage (RSCU) for the mitochondrial protein-coding genes of 


P. lineolatus. 

Codon Aminoacids Number Ratio of Codon RSCU Codon Aminoacids Number Ratio of Codon RSCU 
UUU Phe 28 0.008 0.78 UAU Tyr 69 0.019 0.91 
UUC 44 0.012 1.22 UAC 83 0.023 1.09 
UUA Leu2 40 0.011 0.71 CAU His 137. 0.038 1.01 
UUG 35 0.01 0.62 CAC 135 0.038 0.99 
CUU Leul 59 0.016 1.04 CAA Gln 121 0.034 1.38 
CUC 63 0.018 1.11 CAG 54 0.015 0.62 
CUA 73 0.02 1.29 AAU Asn 149 0.042 0.97 
CUG 70 0.02 1.24 AAC 158 0.044 1.03 
AUU Ile 67 0.019 0.77 AAA Lys 122 0.034 1.46 
AUC 107 0.03 1.23 AAG 45 0.013 0.54 
AUA Met 70 0.02 1.36 GAU Asp 32 0.009 0.85 
AUG 33 0.009 0.64 GAC 43 0.012 1.15 
GUU Val 20 0.006 0.87 GAA Glu 38 0.011 1.21 
GUC 23 0.006 1 GAG 25 0.007 0.79 
GUA 20 0.006 0.87 UGU Cys 31 0.009 1.05 
GUG 29 0.008 1.26 UGC 28 0.008 0.95 
UCU Ser2 70 0.02 1.28 UGA Trp 26 0.007 0.8 
UCC 51 0.014 0.93 UGG 39 0.011 12 
UCA 76 0.021 1.39 CGU Arg 22 0.006 0.75 
UCG 28 0.008 0.51 CGC 41 0.011 1.39 
SU Pro 120 0.033 1.43 CGA 26 0.007 0.88 
CCC 95 0.026 1.13 CGG 29 0.008 0.98 
CCA 87 0.024 1.04 AGU Serl 30 0.008 0.55 
CCG 33 0.009 0.39 AGC 74 0.021 1,35 
ACU Thr 141 0.039 1.39 GGU Gly p 0.008 0.74 
ACC 92 0.026 0.9 GGC 46 0.013 1.27 
ACA 137 0.038 T35 GGA 23 0.006 0.63 
ACG 37 0.01 0.36 GGG 49 0.014 1.35 
GCU Ala 41 0.011 1.18 UAA d 0 0 0 
GCC 47 0.013 1.35 UAG ld 0 0 0 
GCA 34 0.009 0.98 AGA x 0 0 0 
GCG 17 0.005 0.49 AGG * 0 0 0 


* Stop (termination) codon excluded. 


Note: Ratio of codon is equal to the number of each codon divided by the total number of codons (3589). 


cannot be found in P. lineolatus, M. fulvius and Bungarus spp. 
A similar finding was reported in previous studies (Dong 
and Kumazawa, 2005; Kumazawa, 2004). At present, the 
mechanism determining the presence or absence of CSB-2, and 
whether it has phylogenetic significance, is still unclear. The 
motif sequence of CSB-1 was highly similar in snake species, 
while the CSB-3 motif was rarely conserved in previous studies 
(Kumazawa et al., 1996; Kumazawa et al., 1998; Dong and 
Kumazawa, 2005; Douglas et al., 2006). Hence, the recognition 
of CSB-3 was very difficult or ambiguous. Comparison of the 
CSB-3 motif of Naja atra (Qian, 2018) and the two CRs of P. 
lineolatus and the other Elapidae snakes (Figure 3; Table S4) 
revealed that the CSB-3 motif is identical in P. lineolatus and 
Naja spp. but slightly different in other Elapidae. Notably, the 
CSB-3 motif of N. atra was located between T'ASI and CSB-1 to 
form a CSB-3- CSB-1- CSB-2 order (Qian, 2018) rather than the 


common CSB-1- CSB-2- CSB-3 structure. This arrangement 
was different than that reported in previous studies of snakes 
(Kumazawa et al., 1998; Dong and Kumazawa, 2005; Douglas 
et al, 2006) and other vertebrates (Sbisa et al., 1997). The TAS 
motifs of the eight species (incorporating TASI and TAS2) were 
almost identical to those reported in other snakes (Kumazawa 
et al., 1998). No tandem repeats were discovered, while a poly 
T sequence (5-TTTTTTT-3) was detected at both CRs within 
Elapoidea. These length variations are only seen in regions of 
simple repeats of a nucleotide or dinucleotide, suggesting that 
replicational slippage (Levinson and Gutman, 1987) plays a role 


in creating length mutations. 


3.5. Ka/Ks ratio The Ka/Ks ratios of the 13 PCGs varied from 
0.027 to 0.680, but were always less than 1, indicating that these 
genes are evolving under purifying selection (Figure 4). Within 
all compared mitogenomes, the highest average of Ka/Ks ratio 


Asian 
Herpetological 


Research 
302 | Vol. 12 
Alanine ‚A Arginine ¿£ Asparagine y Aspartic oh Cysteine „U 
(Ala) SS (Arg) UA (Asn) 3.8 (Asp) Au (Cys) AU 
G-C UU CG ora c-G 
cS ee UG 66 AC c-G 
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Figure 2 Putative secondary structures of the 22 tRNAs identified in the mitogenome of Psammophis lineolatus. The tRNAs are labeled 
with their corresponding three-letter amino acid codes. Dashed line (—) indicates Watson-Crick base pairing and solid dots (+) between 
Gand U pairs mark canonical base pairings appearing in RNA. 
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Figure 4 Non-synonymous (Ka) and synonymous (Ks) substitution ratio of comparisons between Psammophis lineolatus and the other 
three Elapidae species. BM: Bungarus multicinctus; NN: Naja naja; OH: Ophiophagus hannah; PL: P. lineolatus. 


was observed for ATP8, while COX1 had the lowest ratio. 
ATP8 was indicated to be under positive or relaxed selection 
pressure, which has also been reported in other vertebrates 
(Sun et al, 2020). On the other hand, the low Ka/Ks ratio of 
COX1 is indicative of strong purifying selection, which is in 
accordance with results reported in Python snakes (Dubey et 
al, 2012). Importantly, the Ks value of COX2 for P. lineolatus-B. 
multicinctus (PL-BM) was not applicable. This arose because the 
Jukes and Cantor correction cannot be calculated when the 
proportion of differences is equal or higher than 0.75 (Rozas et 
al, 2017). 


3.6. Phylogenetic analyses Both BI and ML approaches 
provided identical and well-supported tree topologies. Thus, 
only the BI tree is presented, which includes PP as well as 
UFBoot from ML (Figure 5). With limited taxon sampling, 
monophyly of the Caenophidia (advanced snakes) was 
recovered with strong support (PP = 1.0; UFBoot = 100), which 
is largely consistent with results of recent studies (Greiner et 
al, 2019; Pyron et al., 2011; Pyron et al, 2014; Reeder et al., 2015; 
Wiens et al, 2012; Zheng and Wiens, 2016). In the Caenophidia, 
Acrochordidae was recovered as the sister taxon, followed 


successively by Xenodermatidae, Viperidae, Homalopsidae, 
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Figure 5 A majority-rule consensus tree inferred from Bayesian inference using MrBayes v3.2, based on the combined data set of RNA 
genes and protein-coding genes. The phylogenetic position of Psammophiidae is highlighted in gray. Numbers beside the nodes indicate 
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Elapoidea and Colurbroidae. These results also align with recent 
studies (Streicher and Wiens, 2016; Zaher et al., 2019). Monophyly 
of both the Elapoidea and Colurbroidae superfamilies was 
recovered with strong support (PP = 1.0; UFBoot = 100). Notably, 
in the sampled lineages, Psammophiidae was clustered into the 
Elapoidea superfamily (PP = 1.0; UFBoot = 100). Nevertheless, 
according to recent phylogenetic studies (Figueroa et al., 2016; 
Kelly et al., 2009; Zaher et al., 2019), which were based on 
analyses of a few mitochondrial and/or nuclear genes, the 
sister taxon of Psammophiidae is unclear. Thus, it is necessary 
to add more taxa and mitogenome sequences to explore the 
phylogenetic relationships of the advanced snakes in general, 
and to resolve the phylogenetic position of Psammophiidae in 
Elapoidea in particular. 
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Appendix 


Table $1 PCR primers for amplification of the complete mitochondrial genome of Psammophis lineolatus. 


Annealing Approximate 
o. Primer name 5'—3' sequence temperature (°C) product length References 
used in the PCR (bp) 


Z 


16S-R2 CCATYAGGATGTCCTGATCC 


ND1-P3r AATCGHGGGTAKGANGCTCG 


Met-R2 GGTATGGGCCCGARMTAGCTT 


ND2-P2r GAKGAGAAWGCTATTAGYTTTCG 


CO1-P3r GGATTGADGADGCDCCTGC 


CO2-Plr GTAGGAAGAYRACYTCTTCTAT 


CO2-P2r TGSGGCATTTCACTAAGGAA 


CO3-P1r CATGGGCTBGGGTCWACWAGATG 


B169R GATCATGTAATTGTTGCTCCG 


ND3-P1r GGRTCRAAGCCGCAYTCRTA 


GLr GTTCTTCTGAGTCGAAATCAG 


ND4-P2r TTTGCTAGRCAGAAWAGTGC 


ND5-P1r ACAACTATKGTRCTRGARTG 


ND5-P2r ATDGTGTCTTTTGARTAGAAMCC 


Cytb-Plr AGGCAGGTTAGTAGTATGGA 


rGlu-1L GAAAAACCRCCGTTGTWATTCAACTA 


Cytb-P3r TTDGAGAARTTTTCTGGGTC 


E1417R GTATATGCCAGGCCAAGAA 


rCONT-4H CTCGKTTTWGGGGTTTGRCGA 


CHr CTGTTCGATCGCATCTACT 
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Table S2 Summary of taxonomic groups used in this study, including GenBank accession numbers 


Family Species Accession number References 


Acrochordidae Acrochordus granulatus AB177879 Dong and Kumazawa (2015) 


Viperidae Causus defilippi GU045452 Castoe et al. (2009) 


Azemiops feae KJ872487 Geng and Yan (Unpubl. data) 


Trimeresurus albolabris KF311102 Song et al. (2015) 


Protobothrops dabieshanensis KF003004 Huang et al. (2014) 


Crotalus horridus HM641837 Hall et al. (2013) 


Gloydius saxatilis KM434236 Xu et al. (2015a) 


Gloydius brevicaudus EU913477 Chen and Zhao (2009) 


Homalopsidae Enhydris plumbea DQ343650 Yan et al. (2008) 


Elapidae Bungarus fasciatus EU579523 Chen and Zhao (2009) 


Micrurus fulvius GU045453 Castoe et al. (2009) 


Naja kaouthia LC431744 Singchat et al. (2019) 


Naja atra EU913475 Chen and Zhao (2009) 


Sibynophis chinensis KF360246 Oh et al. (2014) 


Cyclophiops major KF148620 Sun et al. (2017) 


Lycodon semicarinatum AB008539 Qian et al. (2018) 


Lycodon flavozonatum KR911720 Qian et al. (2018) 


Elaphe poryphyracea GQ181130 Lin et al. (Unpubl. data) 


Oocatochus rufodorsatus KC990020 Li et al. (2014) 


Pantherophis guttatus AM236349 Douglas et al. (2006) 


Elaphe schrenckii KP888955 Liu and Zhao (2015) 


Elaphe davidi KM401547 Xu et al. (2016a) 


Natricidae Hebius vibakari KP684155 Xu et al. (2016b) 


Dipsadidae Oligodon ningshaanensis KJ719252 Wang (2014) 


Thermophis zhaoermii GQ166168 He et al. (2010) 


Leptodeira septentrionalis EU728590 Mucahy and Macey (2009) 


EU728583 Mucahy and Macey (2009) 


Sibon nebulatus 


Hypsiglena slevini EU728584 Mucahy and Macey (2009) 


Continued Table S2 


Family Species Accession number References 
Hypsiglena jani EU728592 Mucahy and Macey (2009) 
Hypsiglena torquata EU728591 Mucahy and Macey (2009) 
H. ochrorhyncha ochrorhyncha EU728578 Mucahy and Macey (2009) 
H. o. klauberi EU728589 Mucahy and Macey (2009) 
H. o. nuchalata EU728581 Mucahy and Macey (2009) 
Hypsiglena sp. EU728580 Mucahy and Macey (2009) 
Hypsiglena unaocularus KJ486458 Mucahy et al. (2014) 
H. chlorophaea chlorophaea EU728593 Mucahy and Macey (2009) 
H. c. catalinae KJ486459 Mucahy et al. (2014) 
H. c. deserticola EU728587 Mucahy and Macey (2009) 
Outgroup 
Aniliidae Anilius scytale FJ755180 Castoe et al. (2009) 
Tropidophiidae Tropidophis haetianus FJ755181 Castoe et al. (2009) 
Boidae Boa constrictor AB177354 Dong and Kumazawa (2015) 
Cylindrophiidae Cylindrophis ruffus AB179619 Dong and Kumazawa (2015) 
Xenopeltidae Xenopeltis unicolor AB179620 Dong and Kumazawa (2015) 
Pythonidae Python regius AB177878 Dong and Kumazawa (2015) 
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Table S3 Best-fit models and partitioning schemes selected by PartitionFinder 2. 


Subset Best model Partition identity 


ND5. 1st, NDAL. 1st, ATP6. 1st, ATP6. 2nd, NDAL. 3rd, ATP6. 3rd, ATPS8. 3rd, tRNA. Leul, COXI. Ist, 
COXI1. 2nd, ATPS8. 2nd, ND4. 3rd, ND3. 3rd COX2. 3rd, COX3. 2nd, COX3. 3rd, tRNA. Lys, Cytb. 1st, 
pi GTR+I+G NDA. 2nd, Cytb. 2nd, tRNA. Ile, tRNA_Ser1, tRNA_Thr, ND1_1st, ND1_2nd, ND1_3rd ND2. 1st, ND2_2nd, 
ND2. 3rd, ND3. 1st, ND6. 1st, ND6. 2nd, ND6. 3rd, 16S rRNA, tRNA. Phe, 12S rRNA, tRNA. Val, tRNA. Ser2, 
tRNA. Gln, COX2. 1st, ND5. 2nd, COX3. 1st, Cytb. 3rd, NDAL. 2nd, COXI. 3rd, ND3. 2nd, NDA. 1st 


ATPS8. 1st, ND5. 3rd, tRNA. Trp, tRNA. Leu2, tRNA. Met, tRNA. Ala, tRNA. Asn, tRNA. Asp, tRNA. Glu, 


E CIRES tRNA_Gly, tRNA. His, tRNA. Arg, tRNA. Pro 


P3 HKY+I+G COX2_2nd 
P4 SYM+G tRNA_Tyr, tRNA. Cys 


Table $4 Structure traits of the mitochondrial DNA complete control region (D-loop, CR) in Psammophis lineolatus and seven snake species of 
the family Elapidae. 


(a) Psammophis lineolatus 


Unit From ee To Length (bp) 5'-3' motif sequence 

CRII 3560 4571 1012 

C-rich 3601 3620 20 CCCCCCCCCTACCCCCCCCC 
TAS1 3676 3685 10 ACATTAATGG 
TAS2 3717 3726 10 ATAATAATTA 

CSB III 3864 3877 14 INN NIKACIRC IVC 

CSB I 4220 4242 23 TTCCTTAATGCTTGGTAGACATA 
CRI 16 145 17 166 1022 

C-rich 16 191 16 210 20 CCCCCCCCCTACCCCCCCCC 
TAS1 16 266 16 275 10 ACATTAATGG 

TAS2 116 307 16 316 10 ATAATAATTA 

CSB III 16 454 16 467 14 AAAISCHCHANSCH 

CSBI 16 813 16 835 23 TTCCTTAATGCTTGGTAGACATA 
(b) Naja naja 

Unit Tom an To Length (bp) 5'-3' motif sequence 

CRII 3587 4614 1028 

C-rich 3610 3627 18 CCCCCCCCCTCCCCCCCC 
TAS1 3699 3709 11 ACATTAATGAT 

CSB III 3886 3899 14 AAATCCTCTATCCT 

CSB I 4241 4264 23 TTCCTCCATGATTATAAGACATAT 
TAS2 4400 4410 11 AATATAATTAT 

CSB II 4537 4550 14 AATCCCCCCGACAC 

CRI 16 185 17 213 1029 

C-rich 16 205 16 222 18 (69996 9G 91i 9 
TAS1 16 294 16 304 11 ACATTAATGAT 

CSB III 16 481 16 494 14 AAATCCTCTATCCT 

CSB I 16 836 16 859 23 TTCCTCCATGATTATAAGACATAT 
TAS2 16 995 17 005 11 AATATAATTAT 

CSB II 17 132 17 145 14 AATCCCCCCGACAC 

(c) Naja kaouthia 

Unit From Position To Length (bp) 5'-3' motif sequence 

CRU 3581 4608 1028 

C-rich 3604 3621 18 CCCCCCCCCTCCCCCCCC 
TAS1 3693 3703 11 ACATTAATGAT 

CSB III 3880 3894 14 AAATCCTCTATCCT 

CSB I 4236 4259 23 TTCCTCCATGATTATAAGACATAT 
TAS2 4395 4405 11 AATATAATTAT 

CSB II 4531 4544 14 VAVASTIS (9 (9 (8 (667 (8/8 

CRI 16 177 17 203 1027 

C-rich 16195 16212 18 (eee iG e 
TAS1 16 284 16 294 11 ACATTAATGAT 

CSB III 16 472 16 485 14 AAATCCTCTATCCT 

CSB I 16 827 16 850 23 TTCCTCCATGATTATAAGACATAT 
TAS2 16 986 16 996 11 AATATAATTAT 


CSB II 17 123 17 136 14 AATCCCCCCGACAC 


Continued Table S4 


(d) Naja atra 
Unit —. rnrn en To Length (bp) 5'-3' motif sequence 
CRII 3587 4615 1029 
C-rich 3610 3628 19 CCCCCCCCCCTCCCCCCCC 
TAS1 3700 3710 11 ACATTAATGAT 

CSB III 3887 3900 14 AAATCCTCTATCCT 

CSB I 4242 4265 223) TTCCTCCATGATTATAAGACATAT 
TAS2 4401 4411 11 AATATAATTAT 

CSB II 4538 4551 14 AATCCCCCCGACAC 

CRI 16 186 17 216 1031 

C-rich 16 206 16 225 20 CECCCEECECECCTIECECEECEE 
TAS1 16 297 16 307 11 ACATTAATGAT 

CSB III 16 484 16 497 14 ANA TCCNCTANECT 

CSB I 16 839 16 862 23 TTCCTCCATGATTATAAGACATAT 
TAS2 16 998 17 008 11 AATATAATTAT 

CSB II 17 135 17 148 14 AATCCCCCCGACAC 

(e) Ophiophagus hannah 

Unit From Position To Length (bp) 5'-3' motif sequence 

CRII 3572 4589 1018 

C-rich 3594 3613 20 CCCCCCCCCTACCCCCCCCC 
TAS1 3686 3696 11 ACATTAATGAT 

TAS2 3726 3736 11 ACTATAATTAT 

CSB III 3873 3886 14 AAACCCTCTATCCT 

CSBI 4228 4250 23 TTCTTCCATGCTCGTTAGACATA 
CSB II 4416 4429 14 ACACCCCCTGACAC 

CRI 16 260 17 267 1008 

C-rich 16 281 16 300 20 (voco eA 
'TASI 16 373 16 383 11 ACATTAATGAT 

IAS 16 413 16 423 11 ACTATAATTAT 

CSB III 16 560 16 573 14 AAACCCTCTATCCT 

CSBI 16915 16 937 23 TTCTTCCATGCTCGTTAGACATA 
CSB II 17 103 17 116 14 ACACCCCCTGACAC 

(f) Micrurus fulvius 

Unit Ew Position To Length (bp) 5'-3' motif sequence 

CRII 3551 4705 1155 

C-rich 3575 3591 17 CCCCCCCCCTCCCCCCC 
TAS1 3664 3674 11 ACATTAATGTT 

CSB III 3851 3864 14 AAACCCTCTATCCT 

CSBI 4205 4227 23 TTCCTTCATGCTTGGTAGACATA 
TAS2 4586 4596 11 ATAAAAATTAT 

CRI 16 312 17 506 1195 

C-rich 16 335 16 352 18 CCCCCCCCCTCCCCCCCC 
TAS1 16 425 16 435 11 ACATTAATGTT 

CSB III 16 612 16 625 14 AAACCCTCTATCCT 

CSB I 16 966 16 988 25 TTCCTTCATGCTTGGTAGACATA 


TAS2 17 374 17 357 11 ATAAAAATTAT 


Continued Table S4 


(g) Bungarus multioinctus 


Pod 
nit €— e  ————— Length (bp) 5'-3' motif sequence 


c 
md 
4 
o 
3 
4 
o 


O 
to 
N 


3576 3597 CCCCCCCCCTACCCCCCCCGCC 


^ 
re 
Q 
E 


o 
Wn 
u 
= 
=) 


3855 3868 


m 
> 


AAACCCTCTATCCT 


TAS2 4345 4355 


fer 
= 


ATAATAAAACT 


C-rich 16 125 16 148 24 CCCCCCCCCCCTACCCCCCCCGCC 

TAS] 16217, 16227 ACATTAATCTA | 
CSB III 16 406 16 419 14 AAACCCTCTATCCT 

CSBI — — 1670 —— 1678 — 23  — TTCCTTCATGCTTGTTAGACATA ———— 
TAS2 16 896 16 906 11 ATAATAAAACT 

(h) Bungarus fasciatus 

Unit Pu Position To Length (bp) 5'-3' motif sequence 


C-rich 3585 3604 CCCCCCCCCTTCCCCCCCCC 


N 
o 


AAACCCTCTATCCT 


Q 
n 
w 
= 
[a] 


3865 3878 


mi 
4 


TAS2 4278 4288 


e 
= 


ATAGTAATTTT 


2 
d. 
a 

2m 


16 160 16 180 


N 
ar 


CCCCCCCCCCTTCCCCCCCCC 


AAACCCTCTATCCT 


O 
e 
w 
= 
E 


16 441 16 454 


m 
E 


TAS2 16 854 16 864 ATAGTAATTTT 


= 
un 


